
def doMatrixSolve(A,B,tol=10**(-7)): #A must be square n by n matrix, B is 1 by n vector
  x_new=[[1. for i in range(len(B[0]))]]
  converged=False
  while not converged:
    x_old=[x_new[0][:]]
    x_new=gaussSiedel(A,B,x_old)
    
    converged=True #succeeds by default
    for i in xrange(len(x_old[0])):
      if abs(2*(x_new[0][i]-x_old[0][i])/(x_new[0][i]+x_old[0][i]))>=tol:
        converged=False
        break
  return x_new

def gaussSiedel(A,B,x):
  n=len(x[0])
  y=[[0. for l in range(n)]]
  for i in range(len(x[0])):
    val=0 #dummy variable
    for j in range(i):
      val+=A[i][j]*y[0][j]
      #print 'a',i,j
    for j in range(i+1,n):
      val+=A[i][j]*x[0][j]
      #print 'b',i,j
    y[0][i]=(B[0][i]-val)/A[i][i]
  #print y
  return y
      


